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ABSTRACT 

The upwind difference schemes of Godunov, Osher, Roe and van Leer are 
able to resolve one-dimensional steady shocks for the Euler equations within 
one or two mesh intervals. Unfortunately, this resolution is lost in two 
dimensions when the shock crosses the computing grid at an oblique angle. To 
correct this problem, we develop a numerical scheme which automatically 
locates the angle at which a shock might be expected to cross the computing 
grid and then constructs separate finite difference formulas for the flux 
components normal and tangential to this direction. 

We present numerical results which Illustrate the ability of this new 
method to resolve steady oblique shocks. 
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1. Introduction 


Beginning with the work of Godunov [2], considerable effort has been 
expended seeking numerical methods which can solve the equations of gas 
dynamics accurately and sharply resolve discontinuities. 

In the case of one-dimensional flows containing steady discontinuities , 
this effort has paid off. Work by van Leer [14,15,16], Roe [12], Osher 
[1,11], Harten [5] and others Indicates that the mechanism of "shock capture" 
is well understood and that it is now possible to design second ordef accurate 
schemes which resolve steady discontinuities within one or two mesh Intervals 
without wiggles. A recent paper by Harten and Hyman [3] indicates that these 
results can be extended to one dimensional flows with moving discontinuities. 

In the case of two dimensional flows, some progress has been made but 
there is room for improvement. It is well known that the one-dimensional 
results cited above can be reproduced in two-dimensional calculations if the 
computing grid is chosen so that only one set of grid lines crosses a 
discontinuity in the flow. This would appear to be the optimum way to account 
for the two-dimensional aspects of the flow but in practice it is very 
difficult to choose an appropriate grid either a priori or adaptively during a 
calculation. We shall not examine this approach here. 

An alternative approach is to use dimensional splitting with a high 
resolution one-dimensional scheme in each split step. This method Ignores the 
two-dimensional orientation of discontinuities but appears to work 
surprisingly well (cf. Yee et al. [18], Woodward and Collela [17]). Since 
this method is also simple to program, it is quite attractive for practical 
application. 

Despite these attractions, we attempt to do better. In particular, in 
this paper, we derive a two-dimensional scheme which can sharply resolve weak 
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shocks which are skew to the fixed computing grid. Splitting methods have 
difficulty dealing with this situation. An outline of the paper follows. 

In Section 2, we briefly discuss the ideas behind the Godunov [2] method 
and show heurlstlcally why it resolves steady shocks so well. This is 

followed by an introduction to Godunov-type methods with particular emphasis 
on the flux vector splitting method of van Leer [16] which we later use. 

In Section 3, we derive the numerical scheme. First we review the 
rotated difference scheme which Jameson [6] applied to the transonic full 
potential equation. Then we explain how rotated differences, properly 
applied, can Improve shock resolution for the Euler equations. Finally, we 
construct normal and tangential fluxes for the Euler equations and show that 
the resulting method is in conservation form and is consistent with the 
conservation laws. 

In section 4, we consider the choice of angle for the rotated 
differences. We present a choice which seems to be appropriate for shock 
resolution with the Euler equations and then we briefly discuss more general 
discontinuities and more general conservation laws. 

Section 5 contains numerical computations using this scheme and 
comparisons with other schemes. 

Section 6 sumarlzes the present work, discusses the results of the 
present work and presents an outline for future work. 


2. Godunov-type Methods and One Dimensional Problems 

In this section we study the Godunov method for one-dimensional systems 
of conservation laws 

u^ + f(u)^ = 0 


( 2 . 1 ) 
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In particular we wish to determine which features of this scheme are 
responsible for its ability to resolve steady discontinuities. This 
discussion follows closely that of Harten, Lax and van Leer [4], so the reader 
is encouraged to consult their paper for additional details. 

The construction of Godunov's scheme is as follows. At discrete time 

levels tj^, n = 0,1,..., the numerical approximation v(x,t^) to the solution 
u(x,tjj) of (2.1) is taken to be a piecewise constant function of x, le. 

v(x,tn) = Vj for xel^ = ( ( j- V 2 )Ax, ( j+ V 2 )Ax) (2.2) 

To calculate the numerical approximation at the next time level 
t^^^ = t^ + At we first solve exactly the initial value problem 

u + f(u) = 0, 

t X 

u(x,t^) = v(x,tj, (2.3) 

for - «> < X < “ and t < t < t + At: and denote its solution by 

u^(x,t). This means that at each discontinuity of v(x,tj^) we must solve a 
local Riemann problem. The solution to the Riemann problem at say the 
Interface between Ij and ® similarity solution which depends only 

on the states Vj and the ratio (x-( j+l/2)Ax)/(t-t^). We denote 

this solution by u{ (x-( j+l/2)Ax)/(t-t^);Vj,v^j^} . Since signals propagate 

with finite velocity, 

u{(x-(j+l/2)Ax)/(t-tJ;v”,v^j^} == Vj 


when 
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(x-( j+l/2)Ax)/(t-t ) < a 

n JL 

and 

u{(x-(J+l/2)Ax)/(t-t^);Vj,v"^^} = 

when 

(x-( j-l/2)Ax)/[t-t^) > a^ (2.4) 


where aj^ and a^ are the smallest and largest signal velocity respectively. 

If we keep [ la^^^lAt) /Ax < 1/2 where l^max^ largest signal 
speed, there will be no Interaction between neighboring Riemann problems and 
Un(x,t) can be written exactly as 


for 


Uj^(x:,t) = u{(x-( j+l/2)Ax)/(t-tJ ;Vj,Vj_j_j^} 


jAx < X < (j+l)Ax, t^ < t < t^^^ 


(2.5) 


To obtain a piecewise constant approximation v(x,t^^^) for the next time 
step, Godunov averages Ujj(x,t^^.]^), i.e. he sets 


= 1/Axf u (x,t ,i)dx 

j j n^ n+1^ 

j 


( 2 . 6 ) 


In terms of the solution to the local Riemann problems, we can rewrite this as 


1 /A 

Vj = 1/Ax_ 


Ax/2 


:/ u(x/At;v’^_, ,v")dx + 1/Ax/ u( x/A t; v”,v”'*'^) dx (2.7) 

0 J ^ J _^x/2 J J 


Since u^^ is an exact solution to the conservation laws (2.1), we can 

rH"! 

evaluate the integral defining v^ in (2.6) by integrating (2.1) over 


Ijx(tn,tn+l) to get 
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I ^ u(x,t )dx + / f(u (( j+ V 2 )Ax,t))dt 

I . I . t 

3 J n 

^n+1 

- / f(u^(( j- V 2 )Ax,t))dt = 0 (2.8) 

*"n 


or 


where 


n+1 




(2.9) 


>1/2" 


( 2 . 10 ) 


This shows that the Godunov method Is in conservation form with numerical flux 
given by 

f(v,m) = f(u(0;v,o))) (2.11) 


Next we examine how this scheme resolves stationary shocks. For the sake 
of clarity we consider only a single conservation law. Qualitatively similar 
results can be derived for systems of conservation laws (see Lax [8]) but the 
details are far more complicated. 

A shock is an exact discontinuous solution to the Rleraann problem of the 


form 


u 

n 


(x,t) 



X < St 
X > St 


( 2 . 12 ) 


where s satifies the Ranklne-Hugonlot jump condition 


s(u^-Ul) = f(u^) - f[u^) 


(2.13) 


and the entropy condition 
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where 


Hl > S > Hr 


(2.14) 



(2.15) 


A steady shock is a shock with s=0. 

If the initial conditions represent a steady shock located on the 
boundary of a cell, ie. 




for 

xel . 
3 

and 

3 < 

0 

\ 

for 

xel , 
J 

and 

3 > 

0 


(2.16) 


then equation (2.6) will be satisfied exactly for all time levels and (2.16) 
is an exact solution to Godunov's method. 

If the initial conditions represent a steady shock located within a cell, 
say Iq, then equation (2.6) will yield 



for 

3 

< 

0 


u 

m 

for 

3 

= 

0 

(2.17) 


for 

3 

> 

0 



where Ujjj, is some intermediate value between u^ and Uj^. 

At time t^^ the solution to the Riemann problem at the interface 
between I_j^ and Iq will consist of a shock moving at the speed 


If the shock is weak 


% " \ 


(2.18) 


f(uj “ f(uR) +^(ur)(u^-Ur) « f(uR) + aR(u^-UR) 


( 2 . 19 ) 
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Substitute this into (2.18) and note that ~ 


^r( u -u^ 1 
m L'' 


( 2 . 20 ) 


Since u^i is between u^ and Uj^, (u^-u^) and will have 

opposite signs and their quotient will be negative. Equation (2.14) with 
s=0, shows that a^ is negative. Thus S^ is positive, the shock moves to 
the right and by (2.10) 


v.l/ - \ 


( 2 . 21 ) 


Similarly., we see that the solution to the Riemann problem at the 
Interface between Iq and consists of a shock moving at the speed 








m 


If this shock is weak 


( 2 . 22 ) 


df (u, 1 

f(uj « f(u^) +_^(u^-u^) « f(u^) + a^(u^-uj (2.23) 


Substituting (2.23) into (2.22) and using ~ ^^"r^ obtain 

, - , K-”J 


(2.24) 


Since a-^ is positive, by equation (2.14) with s=0 and 
(u^-u^) have opposite signs, S^ is negative and the shock moves to the 


left. Thus by (2.10) 
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^1/2= = “r 

Substitution of (2.21) and (2.25) into (2.9) yields 

^0 = ^0 = ^0 ( 2 . 26 ) 

Thus (2.17) is an exact solution to Godunov's method. 

We note in passing that the choice of numerical fluxes 

■ 'K) 

and 

f(Um,Uj^) = f(u^) 

which does not depend on Uj^^, confines the effect of the averaging which 
created Uj^^ to a single cell. This would not be true of a conventional 
finite difference method whose numerical fluxes would depend on all of its 
arguments. Subsequent time steps would propagate the effect of Uj^ through 
the entire domain and thus spread the shock over many cells. 

The construction of solutions of Riemann problems for nonlinear systems 
is a complicated iterative procedure. In addition, equations (2.9) and (2.10) 
show that, although the entire Riemann solution is computed, only its value at 
the cell interface is actually used. For this reason much recent research has 
been devoted to the construction of numerical flux functions which retain the 
shock capturing ability of Godunov's scheme but which are simpler to 
construct. A particularly simple approach to the construction of numerical 
flux functions for the Euler equations is the flux vector splitting of van 
Leer [16], which we now describe. 
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Van Leer splits the flux f(w) into a forward flux f^(w) and a 
backward flux f~(w) which satisfies 

1. f(w) = f^(w) + f~(w) 

2. df^/dw has all eigenvalues > 0; df~/dw has all eigenvalues < 0 

3. f ± (w) is continuous with f'*’(w) = f(w) for Mach numbers 
M > 1 and f (w) = f(w) for Mach numbers M < -1. 

4. The components of f^ and f together must mimic the symmetry of f 
with respect to M (all other state quantitities held constant), that 
is f^(M) = ± fj^(-M) if f^(M) = ± fj^(-M). 

5. df*/dw must be continuous. 

6. df*/dw must have one eigenvalue vanish for |m1 < 1. 

7. f^(M), like f(M), must be a polynomial in M, and of the lowest 

possible degree. 

If the one dimensional fluxes for the Euler equations with ideal gas law 
are considered to be functions of density p, sound speed c and Mach 
number M, the resulting splittings are 

1) mass 

pu = pcM = pc{ V 2 (M+1)}^ - pc{ V 2 (-M+1)}^ (2.27) 

2 ) momentum 

pu^ + P = pc^[m^+1/y) = pc^{ V 2 (M+1)} ^(-^^^+2/y) 

+ pc^{ I /2 (-M+l)}^(- ^^^+2 /y) (2.28) 


3) total energy 


(e+p)u = pc^m( V 2 M^+1/(T~1)) 


(f+ , ) 

^ momentum^ 
2 [y^- 1) (f'*'mass) 


_ Y 


momentum-' 


(f‘ 


(2.29) 


mass 
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The reader is referred to van Leer's paper for a detailed derivation of 
these expressions. In addition, van Leer shows that these expressions satisfy 
the seven conditions stated above and that these conditions are sufficient to 
assure that steady shocks are resolved within two cells. The numerical flux 
based on this splitting is 

f(u,v) = f^(u) + f (v) (2.30) 

In the next section we Incorporate this numerical flux into a two- 
dimensional rotated difference scheme. 


3. Derivation of a Rotated Difference Scheme 

In a supersonic region, the method of Murman and Cole [10] solves the 
transonic potential equation by replacing derivatives in the streamwise 
direction with upwind difference approximations and replacing derivatives in 
the direction normal to the streamlines with central difference 
approximations. This is easy to do when the computing grid is approximately 
aligned with the streamwise and normal directions but is difficult otherwise. 

Jameson [6] overcomes this problem in the following way. He writes the 
transonic potential equation in a coordinate system aligned with and normal to 
the the local streamwise direction. He then expresses the derivatives in the 
streamwise-normal coordinate system in terms of derivatives in the coordinate 
system of his computing grid, making note of which terms come from streamwise 
derivatives and which terms come from normal derivatives. Finally, those 
terms which come from streamwise derivatives are approximated by upwind 
difference formulas and those terms which come from normal derivatives are 
approximated by central difference formulas. This creates a very effective 


method. 
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In the following we derive a method for the Euler equations which is 
based on this rotated difference idea. In particular, we attempt to choose a 
local coordinate system which permits us to apply the one dimensional theory 
of Section 2. This means that our local coordinate directions must be normal 
and tangential to potential shock directions. In Section 4 we show how to 
determine these directions. Here we assume that the directions are known and 
derive difference formulas that are based on these directions. We note in 
passing that shocks are approximately normal to the streamlines in transonic 
flows. Thus our choice of local coordinate system is equivalent to Jameson's 
in this case. 

Consider a local cartesian coordinate system chosen as above. Such a 
coordinate system has coordinates (x',y'). The coordinate system of the 
computing grid, assumed for simplicity to be cartesian, has coordinates 
(x,y)„ Dependent variables measured in the local and global coordinate 
systems are primed and unprimed respectively. This geometry is shown in 
Figure 1. 

Since the Euler equations are Invariant under rotation, they can be 
written immediately in the local coordinate system as 

where 

U' = [p ,pu' ,pv' ,e]'^ 

F' = [pu',pu'^+ p,pu'v' ,(e+p)u']"^ 


G' = [pv' ,pu'v' ,pv'^+p, (e+p)v' ^ 


(3.2) 
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and 


e = p/(y-l) + V 2 p(u^^+v'^) (3.3) 

= u COS0 + V sln9 


V'' = -u sin0 + V COS0 (3.4) 

In the vicinity of a plane steady shock, each of the terms in (3.1) is 
zero independently. We attempt to construct a numerical scheme which has the 
same property. 

If we express the second terra of (3.1) in global coordinates we get 


8F" 

9x" 


cos9 


3F" 

8x 


+ sln0 


9F" 

9y 


(3.5) 


This can be approximated by a finite difference expression 


COS0 

Ax 


[f-( 







+ 


siH§.[F-(u: ,u: .) - F"(r .,V' . )] 

Ay L l,j+l’ i,j^ ^ i,3 


(3.6) 


The conditions that this expression be zero across a steady shock independent 
of Ax and Ay are 


and 


[f'(u 

[f'[u 


i+l,j 


i,j+l 




)-F^( 






)] =0 
)] =0 


(3.7) 
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These relations are satisfied by the Godunov-type methods discussed In Section 
2. For this reason, the F' fluxes should be approximated by Godunov-type 
numerical fluxes. We use the van Leer flux vector splitting but any other 
Godunov-type flux would suffice. 

The third term in equation (3.1) says that G' is constant along lines 
which are parallel to a shock. To approximate this condition on a discrete 
grid we must make an assumption about the behavior of G' between grid 

points. For lack of any better information we assume that G' varies 
linearly between grid points. 

A representative case is shown in Figure 2. We note that 

-tan9 = Ax/il (3.8) 

so 

i/hy = Ax/(-Ay tanG) = (3.9) 

The condition that G' varies linearly between grid points can be written as 


(g 


i+l,J+l ^1+1, j 


)/Ay = (G-- 




or 


- (g 


1+1. j 


- G' 


I slnO 
' Ax 


(G 


1+1. j+1 




COS0 

Ay 


= 0 


(3.10) 


The condition that G' be constant along a line parallel to a shock can be 
written 

G' = G: . (3.11) 

1, J 


If we express the third term of equation (3.1) in global coordinates, we 


get 
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8G" . A 3G" „ 3G" ... 

■T— ^ = -sin6 -r + COS0 (3.12) 

9y 9x 9y 

We allow the numerical fluxes to depend on additional values of U' and 
approximate this expression by 








(3.13) 


We then identify the numerical fluxes with the corresponding flux values in 
expression (3.10) and (3.11). If our linear interpolation assumption is 
correct, the resulting expression (3.13) Is equal to zero. 

We can derive many expressions that are similar to (3.10) and (3.11) and 
linear combinations of these can be used to construct numerical flux functions 
for (3.13). Thus the formulas presented below are not unique and much work is 
needed to determine which formulas are best. 

To make our G' numerical flux look similar to our F' numerical flux, 
we split it up as follows: 




+ 




(3.14) 
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Q is a parameter chosen to assure that the overall method remain stable. 
Current research is devoted to finding an optimal chloce for this parameter 
and much remains to be done. Futher discussion on the choice of Q appears 
in Section 5. Table 1 specifies how and G' depend on the angle 0. 

Until now we have discussed the derivation of numerical fluxes of local 
variables in local coordinate directions. To compute a global solution we 
need numerical fluxes of global variables in global coordinate directions. To 
obtain these we first express the flux of global variables in global 
coordinates in terras of the flux of local variables in local coordinates as 
shown below. 

Case 1. Flux of a scalar variable 


fj^CU) = pu = pu'cos6 - pv'sin0 = F' (U'' )cos0 - G^(U^)sln0 


Case 2. Flux of a vector component 


F 2 (U) = pu + p = p(u'cos0-v'sin0)(u^cos0-v^sln0 ) + p 


2 2 2 2 2 2 
= pu^ cos 0 - 2pu'v'cos0sin0 + pv' sin 0 + p(cos 0+sln 0) 


F 2 <U')cos 0 - G 2 (U^)cos0sln0 - G^(U' )cos0sln0 


+ F2(U")sin 0 


(3.16) 


The numerical flux is computed by replacing the local flux by its 
corresponding local numerical flux. The two cases shown demonstrate that the 
flux of a scalar variable, such as p or e, transforms like a vector or 
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first order cartesian tensor and the flux of a vector component, such as a 
momentum comonent, transforms like a second order cartesian tensor. These 
transformation laws are well known (cf. Segal [13]) and are easily programmed. 

To summarize, if we are given an angle 0, we compute a numerical flux 
as follows: 

(1) Compute velocity components in the local coordinate system. 

(2) Compute normal and tangential numerical flux components in the local 
coordinate system as described above. 

(3) Use the cartesian tensor transformation laws to compute numerical 
flux components of the global variables in the global coordinate 
system. 

We conclude this section with a brief discussion on what it means for a 
scheme to be conservative and consistent and thus show that the scheme 
described here posseses these characteristics. 

A numerical scheme is called conservative if its numerical flux satisfies 
a discrete version of the divergence theorem. For a rectangular domain, this 
means that we must satisfy an identity of the form 


y y r VV2,j'VV2>J ^ ^1 > 3+ Vz , j- V 2 

iii U 


AxAy = 


N 


JJ^N+V2 ,j“^V2 +I /2 ^i, V2^^"' 


(3.17) 


To satisfy this relation, we associate with each cell boundary point of 
the form ^ single numerical flux value ^ 1 + 1/2 j with each 

cell boundary point of the form ^ ® single numerical flux value 

j+1/2' This assures that the expression on the left of (3.17) telescopes 
to give the expression on the right. 
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Since each numerical flux value, depends on an angle 6 , we associate a 
value of 6 with each cell boundary point. The derivation of the tangential 
local flux G' assumed that 6 was constant through a cell. This would 
indicate that we might prefer to associate the angle 0 with cell centers 
rather than cell boundaries. Unfortunately we could not find a way to do this 
and maintain conservation. 

Lax and Wendroff [7] have shown that a numerical scheme for hyperbolic 
conservation laws is consistent with the differential equations if the 
numerical flux reduces to the differential equation flux when all arguments of 
the numerical flux are set equal. That is 

F(u,u,*»*,u) = F(u) (3.18) 

The scheme described here is derived so that (3.18) is satisfied by the 
numerical flux of local variables in local coordinates. Thus (3.18) is 
satisfied by the global numerical flux because both the global differential 
equation flux and the global numerical flux are computed from their respective 
local flux values by the same formulas. 


4 . Cliioice of Direction 

In Section 3 we describe a method designed to resolve shocks by using 
different numerical flux functions in directions normal and tangential to 
shocks. In this section we show how to find these directions. 

It is Important to note that the algorithm we describe does not determine 
whether or not a shock exists. Thus it is not a shock-fitting algorithm. 
Instead the algorithm always assumes that a steady shock exists and computes 
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its normal direction. Our numerical tests Indicate that this approach locates 
the proper direction when it is needed and causes no problems otherwise. 

To determine the direction of a steady oblique shock, we note ,as do 
Gasdynamics texts (cf. Llepmann and Roshko [9]), that a steady oblique shock 
can be studied as a normal shock with a superimposed uniform tangential 
velocity. Thus, if we are given two velocity vectors, we can locate a 
possible shock direction by finding a rotated coordinate frame in which both 
vectors have a common component. 

A simple way to accomplish this was suggested by John Strikwerda. He 
pointed out that since the velocity in the tangential direction does not 
change across the shock, the shock must be normal to the velocity jump. That 
is, given two velocity vectors 


v(l-l,j) = [u(l-l, j),v(i-l, j)] 

and 

v(l,j) = [u(l, j),v(l, j)] 

the shock, if it exists, is normal to the direction of the vector 

6 v(l,J) = [u(l, j)-u(l-l, j),v(l, j)-v(i-l, j)] = [fi u,6 v] (4.1) 

X XX 


Thus the angle 6 used to compute rotated numerical fluxes is 


0 = arctanfd v/6 ul 

^ X X 


(4.2) 


In practice, the velocity components are computed using a first order 
finite difference formula. In the smooth parts of the flow the error in the 
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velocity differences can be of the same order as the velocity differences 
themselves. This can cause the angle 0 to vary wildly in a part of the flow 
field where little is happening. Numerical experiments indicate that this 
degrades the performance of the method. 

To prevent this, we replace the velocity differences in equation (4.2) by 
weighted averages of the velocity differences at a number of points. At 
present we use the averaging 

6^u(l,j) = [A^u(i-l,j) + 4 A^u(l,j) + A^u(l+1, j)] /6 (4.3) 

where 

A^u(i,j) = {[u(i.j-l)-u(i-l,j-D] +4 [u(i,j)-u(l-l,J)] 

+ [u(l,J+l)-u(i-l,J-l)]}/6 

and 

5yU(i.j) = [4yU(i,j-l) + 4 (a^u( 1, j)+A^u(i, j+l)]/6 

where 

AyU(i,j) = { [u(l-l, j)-u(i-l, j-1)] + 4 [u(l, j)-u(i, j-1)] 

+[u(i+l,j)-u(i+l,j-l)]}/6 (4.6) 

Numerical tests Indicate that this procedure locates the proper shock angles 
and provides smooth angle variations. 

The disadvantage of this approach is that it locates steady shocks for 
the Euler equations but does not locate steady contact discontinuities for the 
Euler equations or steady discontinuities for other systems of conservation 
laws. This limitation is also true of the van Leer formulas that we used to 
compute numerical flux components in the normal direction. In the future we 


(4.4) 

(4.5) 
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plan to construct a method which can locate these more general discontinuities 
and to use the method of Roe [12] to resolve them. 

One way that we might do this is to note that for the scalar equation 

Ut + f(u)^ + g(u)y = 0 (4.7) 

the normal flux is continuous across a steady discontinuity but the tangential 
flux is not. Therefore the discontinuity, if it exists, is located in the 
direction of the vector 

6 T = [f(i,j)-f(i-l,j),g(l,j)-g(i-l,j)] =[6 f,6g]. (4.8) 

X X X 


For systems we might take linear combinations of vectors like (4.8). One 
possibility for a system of m equations would be to seek potential 
discontinuities in the direction of the vector 


6F = 


■ m [6fj)(6uj) m (6gj)(6uj) 
J=1 (Su/ ’j=l (6 u.).2 


(4.9) 


Thus far no numerical experiments have been performed using equation (4.9). 

Baines^^^ has proposed another promising approach to the problem of 
locating general discontinuities. He chooses the angle which satisfies a 
discrete approximation to the equation 


T. Baines, 1982, University of Reading, Reading, England, personal 
communication. 
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9G n 

tt-t- = -SlnO-r — + COS0-r — = 0 (^•10) 

8y 9x 9y 

At this time, some details of this procedure need to be worked out and no 
numerical results are available. 


5« Numerical Results 

In this section we present some numerical computations which demonstrate 
the ability of the present first order method to resolve steady, oblique 
shocks. These results are compared with the results of computations using 
state-of-the-art first and second order upwind methods. 

Consider the problem of supersonic flow over a wedge which is illustrated 
In Figure 3. The solution to this problem is a single oblique shock wave. We 
consider only cases where the flow is supersonic everywhere. 

To solve this problem we construct a uniform computing grid aligned with 
the wedge as shown in Figure 4. We specify all variables at the left and top 
boundaries and we extrapolate all variables at the right boundary. These 
boundary conditions are correct at the left and right boundaries but the top 
boundary is overspecified since the normal velocity at this boundary is 
subsonic. Fortunately, for the cases considered here, no signals reach this 
boundary from inside the computational domain and this overspecification 
causes no difficulties. 

At the lower boundary we specify that the velocity normal to the wall be 
zero and use some numerical procedure to specify the remaining variables at 
the wall. We have experimented with a number of numerical boundary condition 
procedures and have found that the following gives the best results, at least 
for the methods and problems considered here. 
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To compute the pressure and density at the wall, we follow Chakravarty 
and Osher [1] and assume that locally, only simple plane waves leave the 
boundary. In this case the relations 


c 


(Y-1). 


const. 


p/p" = 


const. 



(5.1) 


which are satisfied along the characteristics that reach the wall from the 
interior, and the wall boundary condition, v=0, determine the pressure and 
density at the wall. 

Various schemes have been used to determine the tangential velocity at 
the wall. Here we use the condition that, in the steady state the total 
enthalpy along the wall should be constant. This gives a relation 

= y ■ Y p/p u^/2 = const. (5.2) 

along the wall, which permits us to determine u once p and p are known. 

Our initial conditions consist of the shock jump of the exact solution 
oriented at an angle of 45® to the computing grid. 

Figure 5a is a three-dimensional plot of the density over a 10® wedge 
at Mach 2 computed using the first order method of Osher and Solomon [11]. 
Figure 5b is a density profile along the line indicated by an arrow in Figure 
5a. The shock angle is approximately 30® . The 30® shock was very 
difficult for any method to resolve as can be seen by the fact that the shock 
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is spread over many grid points. This is not an acceptable solution. 

To be fair we note that the results shown for the Osher-Solomon method 
are typical of results for first order upwind schemes. All first order upwind 
schemes have difficulty with weak, oblique shocks and the Osher-Solomon method 
has no more difficulty than any other. 

Figures 6a and 6b show the results of computations of the flow over a 
10^ wedge at Mach 2 using the second order upwind method of van Leer 
described in [14] and [15]. The numerical flux for this method is based on 
the flux vector splitting formulas, described in Section 2, applied to a 
second order approximation to the dependent variables. A slope limiting 
function is used to prevent the oscillations which usually occur when shocks 
are computed using second order methods. These results show that the scheme 
is monotonic and that it spreads the shock over approximately five grid points 
in this case. 

Figures 7a and 7b illustrate an ideal situation. These results were 
computed using the rotated finite difference scheme of Section 3 with the 
exact shock angle specified and the parameter Q set to zero. In this case 
the shock is confined to between two and three grid points. In [16] van Leer 
shows that flux vector splitting requires two grid points to resolve steady 
one-dimensional shocks. This leads us to believe that the results shown are 
the best that can be expected for a two dimensional method based on these 
numerical flux formulas. These results also indicate that the assumptions 
used to derive this method are reasonable. 

Unfortunately, when we select the shock angle by the algorithm of Section 
4, the method is unstable when Q==0. We suspect that this is due to the 
special variations of the computed shock angle since these were not taken into 
account when the method was derived. We are presently trying to account for 


this. 
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At this time 
for the parameter 
choice 


is stable. Figures 8a and 8b demonstrate the performance of this choice 
combined with the automatic angle method of Section 4. The shock is confined 
to approximately three grid points. 

The second problem that we study is the regular reflection of a shock 
from a plane wall. The physical situation is shown in Figure 9. 

The computations shown below were made on a 61 x 21 uniform grid. 
Exact values of all variables were specified on the left and top boundaries. 
All variables were extrapolated at the right boundary and the wall boundary 
conditions described above were imposed on the lower boundary. 

The initial conditions imposed consisted of two shocks which had the same 
jumps as the exact solution but were located in the wrong place. 

Figure 10a is a carpet plot of the density for this flow as computed 
using the first order method of Osher and Solomon. Figure 10b is a section of 
this plot taken at y = .25. These pictures show that this method smears the 
shocks so badly that, at y = .25, they cannot be distinguished. 

Figures 11a and 11b are the corresponding density plots for this flow as 
computed using the second order upwind method of van Leer. These results are 
much better than the previous ones. Here we can distinguish two distinct 
shocks. The first one is spread over between four and five mesh points and 
the second one is spread over between six and seven mesh points. 

Finally, in Figures 12a and 12b, we show density plots for this flow 
computed using the first order rotated method. This method spreads the first 


we can draw no specific conclusions as to the proper choice 
Q. Extensive numerical experimentation indicates that the 


^i+V2 ,j 


v'. . + vl , . . . ^ 
= ^ > J ^-+1 > J At 

2 Ax 


(5.3) 
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shock over between two and three mesh points and the second shock over between 
five and six mesh points. There is a slight undershoot in front of the first 
shock. 

At first sight, these results might appear to be less than dramatic but 
we should note that we have been comparing a first order method with one of 
the more sophisticated of second order methods. Under these circumstances, we 
are encouraged by the fact that the first order method consistently resolved 
steady oblique shocks within fewer mesh intervals than the second order 
method. This confirms that our oblique shock model works. With additional 
effort we can make it work better. 


6 . Summary and Discussion 

In this paper we describe a method which determines the orientation of 
possible shock solutions to the Euler equations. In addition we show how tViis 
information can be used to construct a first order upwind method which 
computes solutions with greatly Improved steady shock resolution. Indeed, we 
have shown that our first order method resolves steady shocks within fewer 
mesh Intervals than a sophisticated second order upwind method. 

This work has shown that the shock resolving ability of numerical methods 
for hyperbolic equations can be Improved if the methods take the orientation 
of possible shocks into account. Although the results thus far have been 
encouraging, we wish to improve the present method. Before this can be done 
some questions need to be answered. In particular, we wish to understand how 
the variation in angle affects the performance of the method and how we should 
select the parameter Q in the tangential flux. At present this parameter is 


selected in an ad hoc manner. 


It is also Important that we construct a 
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careful stability analysis for this method and that we study appropriate 
choices for boundary conditions. It might also be helpful to study other 
choices of upwind formulas for the normal flux. 

On a more practical level, we realize that first order methods are 
usually not accurate enough in regions of smoothly varying flow to be useful 
in applications. Therefore we are presently developing a second order 
accurate version of our method. This work will be described in a future 
paper. 
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Figure Captions 


1. Geometry of local and global coordinate systems. 

2. Construction of tangential flux. 

3. Oblique shock problem in physical domain. 

4. Oblique shock problem in computational domain. 

5a. Three dimensional plot of density for oblique shock problem, 

computed using the method of Osher and Solomon. 

5b. Density profile six grid points from the wall for the oblique 

shock problem, computed using the method of Osher and Solomon. 

6a. Three dimensional plot of density for oblique shock problem, 

computed using the second order method of van Leer. 

6b. Density profile six grid points from the wall for the oblique 

shock problem, computed using the second order method of van 
Leer. 

7a. Three dimensional plot of density for oblique shock problem, 

computed using rotatlonally biased differences with exact shock 


angle. 
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Figure 


7b. Density profile six grid points from the wall, computed using 
rotationally biased differences with exact shock angle. 

8a. Three dimensional plot of density for oblique shock problem 

computed using rotationally biased differences with automatic 
angle algorithm. 

8b. Density profile six grid points from the wall computed using 

rotationally biased differences with automatic angle algorithm 

9. Shock reflection problem. 

10a. Three dimensional plot of density for shock reflection problem 
computed using the method of Osher and Solomon. 

10b. Density profile at y = *25 for shock reflection problem 

computed using the method of Osher and Solomon. 

lla. Three dimensional plot of density for shock reflection problem 
computed using the second order method of van Leer. 

llb. Density profile at y = .25 for shock reflection problem 

computed using the second order method of van Leer. 

12a. Three dimensional plot of density for shock reflection problem 
computed using rotationally biased differences with automatic 


angle algorithm. 
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Figure 12b. Density profile at y = .25 for shock reflection problem 

computed using rotationally biased differences with automatic 
angle algorithm. 
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Figure 5b. 
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Figure 6b. 
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Figure 7a. 
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Figure 7b. 
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Figure 8b. 
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